cd  "E:\Replication Data for Media Freedom and Public Trust in Science" 

use Data_Gu_Main, clear

**# 前期准备
***生成控制变量的宏
*个人层面： female  emp  income 
global cov ( i.female i.emp i.income )
*国家层面： loggdppc logpop  log_nobel_pc log_qs500_pc e_peaveduc e_pelifeex
global agg ( c.loggdppc c.logpop c.log_nobel_pc c.log_qs500_pc c.e_pelifeex )
**#所有控制变量
global control ///
edu female emp income loggdppc logpop log_nobel_pc log_qs500_pc e_pelifeex 

*Alternative measures of media freedom exposure
global Mexp_other ///
Mexp_DSP_14after_std Mexp_scoreRSF_14after_std Mexp_conFOTN_14after_std Mexp_FOTPscore_14after_std Mexp_GMF_14after_std

***画图设置
set scheme s2color
grstyle init
grstyle color background white
grstyle set plain,horizontal grid
grstyle set color Set1
grstyle color ci_area gs12%40

**#Summary Statistics
**#Table A.1: Summary Statistics for Survey Data
global sum ///
y2 edu female emp income birthyear loggdppc logpop log_nobel_pc log_qs500_pc e_pelifeex ///
Mexp_conVDem_14after Mexp_binVDem_14after Mexp_DSP_14after Mexp_scoreRSF_14after  Mexp_conFOTN_14after Mexp_FOTPscore_14after Mexp_GMF_14after ///
trust_govt trust_journalists science_literacy 

eststo summ: estpost summarize  $sum 

esttab summ using TA1.rtf,  ///
nomtitle nonumber noobs replace label varwidth(35) modelwidth(9) ///
cells("mean(label(Mean) fmt(a2)) sd(label(SD) fmt(a2)) min(label(Min)) max(label(Max)) count(label(Obs) fmt(0))")  ///
title(\b Summary Statistics for Survey Data)  ///
varlabel( Mexp_conVDem_14after "Continuous Media Freedom Exposure (V-Dem), age 14 to present" Mexp_binVDem_14after "Binary Media Freedom Exposure (V-Dem), age 14 to present" Mexp_DSP_14after "Digital Media Freedom Exposure (V-Dem), age 14 to present" Mexp_scoreRSF_14after "Press Freedom Exposure (RSF), age 14 to present" Mexp_conFOTNscore_14after "Internet Freedom Exposure (FOTN), age 14 to present" Mexp_FOTPscore_14after "Press Freedom Exposure (FOTP), age 14 to present" Mexp_GMF_14after "Media Freedom Exposure (GMF), age 14 to present"  ///
trust_govt "Trust in government" trust_journalists "Trust in media" science_literacy "Science Literacy" ///
birthyear "Birth year" female "Female" emp "Employment status" income "Income" edu "Education level" loggdppc "Log GDP per capita" logpop "Log population" log_nobel_pc "Log number of Nobel Prize winners per capita" log_qs500_pc "Log number of QS 500 universities per capita" e_pelifeex "Life expectancy" ) 

**#Additional Results
**#Table A2.Multilevel Model: Media Freedom Exposure 
eststo clear

local var Mexp_DSP_14after
eststo h1: mixed y2 edu female emp income loggdppc logpop log_nobel_pc log_qs500_pc e_pelifeex `var'_std || country_new: , mle
estat icc
estadd scalar ICC = r(icc2)*100
qui tab country if `var'_std != .
estadd scalar N2 = r(r)

local var Mexp_scoreRSF_14after
eststo h2: mixed y2 edu female emp income loggdppc logpop log_nobel_pc log_qs500_pc e_pelifeex `var'_std || country:, mle
estat icc
estadd scalar ICC = r(icc2)*100
qui tab country if `var'_std != .
estadd scalar N2 = r(r)

local var Mexp_conFOTN_14after 
eststo h3: mixed y2 edu female emp income loggdppc logpop log_nobel_pc log_qs500_pc e_pelifeex `var'_std || country:, mle
estat icc
estadd scalar ICC = r(icc2)*100
qui tab country if `var'_std != .
estadd scalar N2 = r(r)

local var Mexp_FOTPscore_14after
eststo h4: mixed y2 edu female emp income loggdppc logpop log_nobel_pc log_qs500_pc e_pelifeex `var'_std || country:, mle
estat icc
estadd scalar ICC = r(icc2)*100
qui tab country if `var'_std != .
estadd scalar N2 = r(r)

local var Mexp_GMF_14after
eststo h5: mixed y2 edu female emp income loggdppc logpop log_nobel_pc log_qs500_pc e_pelifeex `var'_std || country:, mle
estat icc
estadd scalar ICC = r(icc2)*100
qui tab country if `var'_std != .
estadd scalar N2 = r(r)

esttab h* using TA2.rtf , b(3) se(3) nonote nobaselevel noomitted  eqlabel(none) label  replace  star(+ 0.1 * 0.05 ** 0.01 *** 0.01)  /// 
refcat( edu "{\i Individual-Level Variables}" loggdppc "{\i Country-Level Variables}" , nolabel) /// 
keep( _cons $control $Mexp_other ) ///
order(_cons) ///
varlabel( _cons "Constant" female "Female" emp "Employment status" income "Income" edu "Education level" loggdppc "Log GDP per capita" logpop "Log population" log_nobel_pc "Log number of Nobel Prize winners per capita" log_qs500_pc "Log number of QS 500 universities per capita" e_pelifeex "Life expectancy" Mexp_DSP_14after_std "Media freedom exposure") ///
mgroup("{\i DV = Trust in Science (IRT)}") ///
mtitles("{\i Digital Media Freedom Exposure (DSP)}" "{\i Press Freedom Exposure (RSF)}" "{\i Internet Freedom Exposure (FOTN)}" "{\i Press Freedom Exposure (FOTP)}" "{\i Media Freedom Exposure (GMF)}") ///
title({\b Multilevel Model}) ///
stats( ICC N2 N, labels( "ICC (%)" "Country-Level Observations" "Individual-Level Observations") fmt( 2 0 0 ) ) ///
note("{\i \b Notes:} + {\i p} < 0.1, * {\i p} < 0.05, ** {\i p} < 0.01,*** {\i p} < 0.001 (two-tailed test)." )

**#Table A3.Multilevel Model: Media Freedom Exposure and Education Level
eststo clear
local var Mexp_conVDem_14after
eststo h1: mixed y2 female emp income loggdppc logpop log_nobel_pc log_qs500_pc e_pelifeex c.`var'_std##c.edu || country:, mle
estat icc
estadd scalar ICC = r(icc2)*100
qui tab country if `var' != .
estadd scalar N2 = r(r)

local var Mexp_binVDem_14after
eststo h2: mixed y2 female emp income loggdppc logpop log_nobel_pc log_qs500_pc e_pelifeex c.`var'_std##c.edu || country:, mle
estat icc
estadd scalar ICC = r(icc2)*100
qui tab country if `var' != .
estadd scalar N2 = r(r)

local var Mexp_DSP_14after
eststo h3: mixed y2 female emp income loggdppc logpop log_nobel_pc log_qs500_pc e_pelifeex c.`var'_std##c.edu || country:, mle
estat icc
estadd scalar ICC = r(icc2)*100
qui tab country if `var' != .
estadd scalar N2 = r(r)

local var Mexp_conFOTN_14after 
eststo h4: mixed y2 female emp income loggdppc logpop log_nobel_pc log_qs500_pc e_pelifeex c.`var'_std##c.edu || country:, mle
estat icc
estadd scalar ICC = r(icc2)*100
qui tab country if `var' != .
estadd scalar N2 = r(r)

local var Mexp_FOTPscore_14after
eststo h5: mixed y2 female emp income loggdppc logpop log_nobel_pc log_qs500_pc e_pelifeex c.`var'_std##c.edu || country:, mle
estat icc
estadd scalar ICC = r(icc2)*100
qui tab country if `var' != .
estadd scalar N2 = r(r)

local var Mexp_GMF_14after
eststo h6: mixed y2 female emp income loggdppc logpop log_nobel_pc log_qs500_pc e_pelifeex c.`var'_std##c.edu || country:, mle
estat icc
estadd scalar ICC = r(icc2)*100
qui tab country if `var' != .
estadd scalar N2 = r(r)

esttab h* using TA2.rtf , b(3) se(3) nonote nobaselevel noomitted  eqlabel(none) label  replace  star(+ 0.1 * 0.05 ** 0.01 *** 0.01)  /// 
refcat(edu "{\i Individual-Level Variables}" loggdppc "{\i Country-Level Variables}" , nolabel) /// 
keep( _cons $control Mexp_conVDem_14after_std Mexp_binVDem_14after_std Mexp_DSP_14after_std Mexp_conFOTN_14after_std Mexp_FOTPscore_14after_std Mexp_GMF_14after_std c.Mexp_conVDem_14after_std#c.edu c.Mexp_binVDem_14after_std#c.edu c.Mexp_DSP_14after_std#c.edu c.Mexp_conFOTN_14after_std#c.edu c.Mexp_FOTPscore_14after_std#c.edu c.Mexp_GMF_14after_std#c.edu ) ///
order(_cons edu) ///
varlabel( _cons "Constant" female "Female" emp "Employment status" income "Income" edu "Education level" loggdppc "Log GDP per capita" logpop "Log population" log_nobel_pc "Log number of Nobel Prize winners per capita" log_qs500_pc "Log number of QS 500 universities per capita" e_pelifeex "Life expectancy" Mexp_conVDem_14after_std "Media freedom exposure" c.Mexp_conVDem_14after_std#c.edu "Media freedom exposure × Education level") ///
mgroup("{\i DV = Trust in Science (IRT)}") ///
mtitles("{\i Continuous Media Freedom Exposure (V-Dem)}" "{\i Binary Media Freedom Exposure (V-Dem)}" "{\i Digital Media Freedom Exposure (DSP)}" "{\i Internet Freedom Exposure (FOTN)}" "{\i Press Freedom Exposure (FOTP)}" "{\i Media Freedom Exposure (GMF)}") ///
title({\b Multilevel Model: Lifetime Exposure to Media Freedom Moderated by Education Level }) ///
stats( ICC N2 N, labels( "ICC (%)" "Country-Level Observations" "Individual-Level Observations") fmt( 2 0 0 ) ) ///
note("{\i \b Notes:} + {\i p} < 0.1, * {\i p} < 0.05, ** {\i p} < 0.01,*** {\i p} < 0.001 (two-tailed test)." )

**#Figure A1. Marginal Effect of Media Freedom Exposure as Moderated by Education Level
local var Mexp_DSP_14after_std
local title Digital Media Freedom Exposure (DSP)
local num h2
mixed y2 female emp income loggdppc logpop log_nobel_pc log_qs500_pc e_pelifeex c.`var'##c.edu || country:, mle
margins, dydx(`var') at(edu=(1 2 3))
marginsplot,recast(line) recastci(rarea) ciopts(fcolor(%30) lcolor(%1)) nolabels ///
title("`title'",size(3.5)) ytitle("") xtitle(" ")  ///
xsize(6) ysize(4) ylabel(-0.25(0.05)0.05, labsize(small))
graph save F1_`num'.gph, replace  

local var Mexp_conFOTN_14after_std
local title Internet Freedom Exposure (FOTN)
local num h3
mixed y2 female emp income loggdppc logpop log_nobel_pc log_qs500_pc e_pelifeex c.`var'##c.edu || country:, mle
margins, dydx(`var') at(edu=(1 2 3))
marginsplot,recast(line) recastci(rarea) ciopts(fcolor(%30) lcolor(%1)) nolabels ///
title("`title'",size(3.5)) ytitle("") xtitle(" ")  ///
xsize(6) ysize(4) ylabel(-0.25(0.05)0.05, labsize(small))
graph save F1_`num'.gph, replace 

local var Mexp_FOTPscore_14after_std
local title Press Freedom Exposure (FOTP)
local num h4
mixed y2 female emp income loggdppc logpop log_nobel_pc log_qs500_pc e_pelifeex c.`var'##c.edu || country:, mle
margins, dydx(`var') at(edu=(1 2 3))
marginsplot,recast(line) recastci(rarea) ciopts(fcolor(%30) lcolor(%1)) nolabels ///
title("`title'",size(3.5)) ytitle("") xtitle(" ")  ///
xsize(6) ysize(4) ylabel(-0.25(0.05)0.05, labsize(small))
graph save F1_`num'.gph, replace 

local var Mexp_GMF_14after_std
local title Media Freedom Exposure (GMF)
local num h5
mixed y2 female emp income loggdppc logpop log_nobel_pc log_qs500_pc e_pelifeex c.`var'##c.edu || country:, mle
margins, dydx(`var') at(edu=(1 2 3))
marginsplot,recast(line) recastci(rarea) ciopts(fcolor(%30) lcolor(%1)) nolabels ///
title("`title'",size(3.5)) ytitle("") xtitle(" ")  ///
xsize(6) ysize(4) ylabel(-0.25(0.05)0.05, labsize(small))
graph save F1_`num'.gph, replace 

graph combine F1_h2.gph F1_h3.gph F1_h4.gph F1_h5.gph , ///
cols(2) xsize(15) ysize(12) 

graph export "FA1.png", replace

**# Table A4. Fixed Effects Model: Media Freedom Exposure
eststo clear

local var Mexp_conVDem_14after_std
eststo m1: reghdfe y2     ///
`var'    $agg##i.edu $cov##c.`var'   ,a(country i.birthyear) cluster(country)  keepsing
local r2=e(r2_a)
estadd scalar r2a=`r2'

local var Mexp_binVDem_14after_std
eststo m2: reghdfe y2     ///
`var'    $agg##i.edu $cov##c.`var'   ,a(country i.birthyear) cluster(country)  keepsing
local r2=e(r2_a)
estadd scalar r2a=`r2'

local var Mexp_DSP_14after_std
eststo m3: reghdfe y2     ///
`var'    $agg##i.edu $cov##c.`var'   ,a(country i.birthyear) cluster(country)  keepsing
local r2=e(r2_a)
estadd scalar r2a=`r2'

local var Mexp_scoreRSF_14after_std
eststo m4: reghdfe y2     ///
`var'    $agg##i.edu $cov##c.`var'   ,a(country i.birthyear) cluster(country)  keepsing
local r2=e(r2_a)
estadd scalar r2a=`r2'

local var Mexp_conFOTN_14after_std
eststo m5: reghdfe y2     ///
`var'    $agg##i.edu $cov##c.`var'   ,a(country i.birthyear) cluster(country)  keepsing
local r2=e(r2_a)
estadd scalar r2a=`r2'

local var Mexp_FOTPscore_14after_std
eststo m6:reghdfe y2     ///
`var'    $agg##i.edu $cov##c.`var'   ,a(country i.birthyear) cluster(country)  keepsing
local r2=e(r2_a)
estadd scalar r2a=`r2'

local var Mexp_GMF_14after_std
eststo m7: reghdfe y2     ///
`var'    $agg##i.edu $cov##c.`var'   ,a(country i.birthyear) cluster(country)  keepsing
local r2=e(r2_a)
estadd scalar r2a=`r2'

estadd local ctrl2 "√" , replace: m*
estadd local ctrl1 "√" , replace: m*
estadd local rfe "√" , replace: m*
estadd local bfe "√" , replace: m*

esttab m* using TA4.rtf , b(3) se(3) nonote nobaselevel  unstack noomitted  eqlabel(none) label  replace  star(+ 0.1 * 0.05 ** 0.01 *** 0.01) /// 
varlabel(Mexp_conVDem_14after_std "Media Freedom Exposure") ///
keep( Mexp_conVDem_14after_std Mexp_binVDem_14after_std $Mexp_other ) ///
mgroup("{\i DV = Trust in Science (IRT)}", pattern(1 0 0 0 0 0 0) ) ///
mtitles( "{\i Continuous Media Freedom Exposure (V-Dem)}" "{\i Binary Media Freedom Exposure (V-Dem)}" "{\i Digital Media Freedom Exposure (DSP)}"  "{\i Press Freedom Exposure (RSF)}" "{\i Internet Freedom Exposure (FOTN)}" "{\i Press Freedom Exposure (FOTP)}" "{\i Media Freedom Exposure (GMF)}") ///
stats(rfe bfe ctrl1 ctrl2 r2a N, labels("Country FE" "Birth year FE" "Individual-level controls × Media freedom" "Country-level controls × Education level" "Adjusted R²"  "Observations") fmt( 0 0 0 0 3 0 ) ) ///
title({\b Fixed Effects Model: Media Freedom Exposure}) ///
note("{\i \b Notes:} + {\i p} < 0.1, * {\i p} < 0.05, ** {\i p} < 0.01,*** {\i p} < 0.001 (two-tailed test)." )


**#Table A5. Fixed Effects Model: Alternative Measures of Media Freedom Exposure and Education Level
eststo clear

local var Mexp_DSP_14after
reghdfe y2     ///
c.`var'##i.edu    $agg##i.edu $cov##c.`var'   ,a(country i.birthyear) cluster(country)  keepsing
local r2=e(r2_a)
eststo m1: expint y2  `var'  beta
estadd local ctrl2 "√"
estadd local ctrl1 "√"
estadd local rfe "√"
estadd local bfe "√"
estadd scalar r2a=`r2'
lincom primary-secondary
lincom primary-tertiary
lincom secondary-tertiary

local var Mexp_scoreRSF_14after
reghdfe y2     ///
c.`var'##i.edu    $agg##i.edu $cov##c.`var'   ,a(country i.birthyear) cluster(country)  keepsing
local r2=e(r2_a)
eststo m2: expint y2  `var'  beta
estadd local ctrl2 "√"
estadd local ctrl1 "√"
estadd local rfe "√"
estadd local bfe "√"
estadd scalar r2a=`r2'
lincom primary-secondary
lincom primary-tertiary
lincom secondary-tertiary

local var Mexp_conFOTN_14after
reghdfe y2     ///
c.`var'##i.edu    $agg##i.edu $cov##c.`var'   ,a(country i.birthyear) cluster(country)  keepsing
local r2=e(r2_a)
eststo m3: expint y2  `var'  beta
estadd local ctrl2 "√"
estadd local ctrl1 "√"
estadd local rfe "√"
estadd local bfe "√"
estadd scalar r2a=`r2'
lincom primary-secondary
lincom primary-tertiary
lincom secondary-tertiary

local var Mexp_FOTPscore_14after
reghdfe y2     ///
c.`var'##i.edu    $agg##i.edu $cov##c.`var'   ,a(country i.birthyear) cluster(country)  keepsing
local r2=e(r2_a)
eststo m4: expint y2  `var'  beta
estadd local ctrl2 "√"
estadd local ctrl1 "√"
estadd local rfe "√"
estadd local bfe "√"
estadd scalar r2a=`r2'
lincom primary-secondary
lincom primary-tertiary
lincom secondary-tertiary

local var Mexp_GMF_14after
reghdfe y2     ///
c.`var'##i.edu    $agg##i.edu $cov##c.`var'   ,a(country i.birthyear) cluster(country)  keepsing
local r2=e(r2_a)
eststo m5: expint y2  `var'  beta
estadd local ctrl2 "√"
estadd local ctrl1 "√"
estadd local rfe "√"
estadd local bfe "√"
estadd scalar r2a=`r2'
lincom primary-secondary
lincom primary-tertiary
lincom secondary-tertiary

esttab m* using TA5.rtf , b(3) se(3) nonote nobaselevel  unstack noomitted  eqlabel(none) label  replace  star(+ 0.1 * 0.05 ** 0.01 *** 0.01) refcat(primary "{\i Effect of cumulative exposure to media freedom on}" , nolabel) /// 
varlabel(primary "Primary education or below" secondary "Secondary education" tertiary "College education or above") ///
mgroup("{\i DV = Trust in Science (IRT)}", pattern(1 0 0 0 0 0 0) ) ///
mtitles( "{\i Digital Media Freedom Exposure (DSP)}"  "{\i Press Freedom Exposure (RSF)}" "{\i Internet Freedom Exposure (FOTN)}" "{\i Press Freedom Exposure (FOTP)}" "{\i Media Freedom Exposure (GMF)}") ///
stats(rfe bfe ctrl1 ctrl2 r2a N, labels("Country FE" "Birth year FE" "Individual-level controls × Media freedom" "Country-level controls × Education level" "Adjusted R²"  "Observations") fmt( 0 0 0 0 3 0 ) ) ///
title({\b Fixed Effects Model: Alternative Measures of Media Freedom Exposure}) ///
note("{\i \b Notes:} + {\i p} < 0.1, * {\i p} < 0.05, ** {\i p} < 0.01,*** {\i p} < 0.001 (two-tailed test)." )

**#Table A6. Fixed Effects Model: Comparison of Estimates for Different Education Levels
eststo clear

local var Mexp_conVDem_14after
reghdfe y2     ///
c.`var'##i.edu    $agg##i.edu $cov##c.`var'   ,a(country i.birthyear) cluster(country)  keepsing
local r2=e(r2_a)
eststo m1: expint y2  `var'  beta
estadd local ctrl2 "√"
estadd local ctrl1 "√"
estadd local rfe "√"
estadd local bfe "√"
estadd scalar r2a=`r2'
lincom primary-secondary
lincom primary-tertiary
lincom secondary-tertiary

local var Mexp_binVDem_14after
reghdfe y2     ///
c.`var'##i.edu    $agg##i.edu $cov##c.`var'   ,a(country i.birthyear) cluster(country)  keepsing
local r2=e(r2_a)
eststo m2: expint y2  `var'  beta
estadd local ctrl2 "√"
estadd local ctrl1 "√"
estadd local rfe "√"
estadd local bfe "√"
estadd scalar r2a=`r2'
lincom primary-secondary
lincom primary-tertiary
lincom secondary-tertiary

local var Mexp_DSP_14after
reghdfe y2     ///
c.`var'##i.edu    $agg##i.edu $cov##c.`var'   ,a(country i.birthyear) cluster(country)  keepsing
local r2=e(r2_a)
eststo m3: expint y2  `var'  beta
estadd local ctrl2 "√"
estadd local ctrl1 "√"
estadd local rfe "√"
estadd local bfe "√"
estadd scalar r2a=`r2'
lincom primary-secondary
lincom primary-tertiary
lincom secondary-tertiary

local var Mexp_scoreRSF_14after
reghdfe y2     ///
c.`var'##i.edu    $agg##i.edu $cov##c.`var'   ,a(country i.birthyear) cluster(country)  keepsing
local r2=e(r2_a)
eststo m4: expint y2  `var'  beta
estadd local ctrl2 "√"
estadd local ctrl1 "√"
estadd local rfe "√"
estadd local bfe "√"
estadd scalar r2a=`r2'
lincom primary-secondary
lincom primary-tertiary
lincom secondary-tertiary

local var Mexp_conFOTN_14after
reghdfe y2     ///
c.`var'##i.edu    $agg##i.edu $cov##c.`var'   ,a(country i.birthyear) cluster(country)  keepsing
local r2=e(r2_a)
eststo m5: expint y2  `var'  beta
estadd local ctrl2 "√"
estadd local ctrl1 "√"
estadd local rfe "√"
estadd local bfe "√"
estadd scalar r2a=`r2'
lincom primary-secondary
lincom primary-tertiary
lincom secondary-tertiary

local var Mexp_FOTPscore_14after
reghdfe y2     ///
c.`var'##i.edu    $agg##i.edu $cov##c.`var'   ,a(country i.birthyear) cluster(country)  keepsing
local r2=e(r2_a)
eststo m6: expint y2  `var'  beta
estadd local ctrl2 "√"
estadd local ctrl1 "√"
estadd local rfe "√"
estadd local bfe "√"
estadd scalar r2a=`r2'
lincom primary-secondary
lincom primary-tertiary
lincom secondary-tertiary

local var Mexp_GMF_14after
reghdfe y2     ///
c.`var'##i.edu    $agg##i.edu $cov##c.`var'   ,a(country i.birthyear) cluster(country)  keepsing
local r2=e(r2_a)
eststo m7: expint y2  `var'  beta
estadd local ctrl2 "√"
estadd local ctrl1 "√"
estadd local rfe "√"
estadd local bfe "√"
estadd scalar r2a=`r2'
lincom primary-secondary
lincom primary-tertiary
lincom secondary-tertiary

estimate restore m1
eststo mx1: nlcom (p_s:(_b[primary]-_b[secondary])) ///
				(s_t:(_b[secondary]-_b[tertiary])) ///
				(p_t:(_b[primary]-_b[tertiary])), post
				
estimate restore m2
eststo mx2: nlcom (p_s:(_b[primary]-_b[secondary])) ///
				(s_t:(_b[secondary]-_b[tertiary])) ///
				(p_t:(_b[primary]-_b[tertiary])), post

estimate restore m3
eststo mx3: nlcom (p_s:(_b[primary]-_b[secondary])) ///
				(s_t:(_b[secondary]-_b[tertiary])) ///
				(p_t:(_b[primary]-_b[tertiary])), post

estimate restore m4
eststo mx4: nlcom (p_s:(_b[primary]-_b[secondary])) ///
				(s_t:(_b[secondary]-_b[tertiary])) ///
				(p_t:(_b[primary]-_b[tertiary])), post

estimate restore m5
eststo mx5: nlcom (p_s:(_b[primary]-_b[secondary])) ///
				(s_t:(_b[secondary]-_b[tertiary])) ///
				(p_t:(_b[primary]-_b[tertiary])), post
				
estimate restore m6
eststo mx6: nlcom (p_s:(_b[primary]-_b[secondary])) ///
				(s_t:(_b[secondary]-_b[tertiary])) ///
				(p_t:(_b[primary]-_b[tertiary])), post				

estimate restore m7
eststo mx7: nlcom (p_s:(_b[primary]-_b[secondary])) ///
				(s_t:(_b[secondary]-_b[tertiary])) ///
				(p_t:(_b[primary]-_b[tertiary])), post				

esttab  mx* using TA6.rtf, b(3) se(3) nonote nobaselevel  unstack noomitted  eqlabel(none) label replace star(+ 0.1 * 0.05 ** 0.01 *** 0.01) ///
varlabel(p_s "Primary - Secondary" s_t "Secondary - College" p_t "Primary - College") ///
mgroup("{\i DV = Trust in Science (IRT)}", pattern(1 0 0 0 0 0 0) ) ///
mtitles( "{\i Continuous Media Freedom Exposure (V-Dem)}" "{\i Binary Media Freedom Exposure (V-Dem)}" "{\i Digital Media Freedom Exposure (DSP)}"  "{\i Press Freedom Exposure (RSF)}" "{\i Internet Freedom Exposure (FOTN)}" "{\i Press Freedom Exposure (FOTP)}" "{\i Media Freedom Exposure (GMF)}") ///
title({\b Comparison of Estimates for Different Education Levels: Fixed Effects Model}) ///
note("{\i \b Notes:} + {\i p} < 0.1, * {\i p} < 0.05, ** {\i p} < 0.01,*** {\i p} < 0.001 (two-tailed test)." )

**#Robustness Checks
**#Table A.8: Robustness for Multilevel Model: Alternative Media Freedom Exposure Measures
eststo clear

local var Mexp_conVDem_11after
eststo h1: mixed y2 edu female emp income loggdppc logpop log_nobel_pc log_qs500_pc e_pelifeex `var'_std || country:, mle
estat icc
estadd scalar ICC = r(icc2)*100
qui tab country if `var'_std != .
estadd scalar N2 = r(r)

local var Mexp_conVDem_6after
eststo h2: mixed y2 edu female emp income loggdppc logpop log_nobel_pc log_qs500_pc e_pelifeex `var'_std || country:, mle
estat icc
estadd scalar ICC = r(icc2)*100
qui tab country if `var'_std != .
estadd scalar N2 = r(r)

local var Mexp_binVDem_11after
eststo h3: mixed y2 edu female emp income loggdppc logpop log_nobel_pc log_qs500_pc e_pelifeex `var'_std || country:, mle
estat icc
estadd scalar ICC = r(icc2)*100
qui tab country if `var'_std != .
estadd scalar N2 = r(r)

local var Mexp_binVDem_6after
eststo h4: mixed y2 edu female emp income loggdppc logpop log_nobel_pc log_qs500_pc e_pelifeex `var'_std || country:, mle
estat icc
estadd scalar ICC = r(icc2)*100
qui tab country if `var'_std != .
estadd scalar N2 = r(r)

esttab h* using TA7.rtf , b(3) se(3) nonote nobaselevel noomitted  eqlabel(none) label  replace  star(+ 0.1 * 0.05 ** 0.01 *** 0.01)  /// 
refcat( edu "{\i Individual-Level Variables}" loggdppc "{\i Country-Level Variables}" , nolabel) /// 
keep( _cons $control Mexp_conVDem_11after_std Mexp_conVDem_6after_std Mexp_binVDem_11after_std Mexp_binVDem_6after_std ) ///
order(_cons) ///
varlabel( _cons "Constant" female "Female" emp "Employment status" income "Income" edu "Education level" loggdppc "Log GDP per capita" logpop "Log population" log_nobel_pc "Log number of Nobel Prize winners per capita" log_qs500_pc "Log number of QS 500 universities per capita" e_pelifeex "Life expectancy" Mexp_conVDem_14after_std "Media freedom exposure") ///
mgroup("{\i DV = Trust in Science (IRT)}") ///
mtitles("{\i Continuous Media Freedom Exposure (V-Dem), age 11 to present}" "{\i Continuous Media Freedom Exposure (V-Dem), age 6 to present}" "{\i Binary Media Freedom Exposure (V-Dem), age 11 to present}" "{\i Binary Media Freedom Exposure (V-Dem), age 6 to present}") ///
title({\b Robustness for Multilevel Model: Alternative Measures of Media Freedom Exposure and Education Level}) ///
stats( ICC N2 N, labels( "ICC (%)" "Country-Level Observations" "Individual-Level Observations") fmt( 2 0 0 ) ) ///
note("{\i \b Notes:} + {\i p} < 0.1, * {\i p} < 0.05, ** {\i p} < 0.01,*** {\i p} < 0.001 (two-tailed test)." )

**#Table A.8: Robustness for Fixed Effects Model: Alternative Media Freedom Exposure Measures
eststo clear

local var Mexp_conVDem_11after
reghdfe y2     ///
c.`var'##i.edu    $agg##i.edu $cov##c.`var'   ,a(country i.birthyear) cluster(country)  keepsing
local r2=e(r2_a)
eststo m1: expint y2  `var'  beta
estadd local ctrl2 "√"
estadd local ctrl1 "√"
estadd local rfe "√"
estadd local bfe "√"
estadd scalar r2a=`r2'
lincom primary-secondary
lincom primary-tertiary
lincom secondary-tertiary

local var Mexp_conVDem_6after
reghdfe y2     ///
c.`var'##i.edu    $agg##i.edu $cov##c.`var'   ,a(country i.birthyear) cluster(country)  keepsing
local r2=e(r2_a)
eststo m2: expint y2  `var'  beta
estadd local ctrl2 "√"
estadd local ctrl1 "√"
estadd local rfe "√"
estadd local bfe "√"
estadd scalar r2a=`r2'
lincom primary-secondary
lincom primary-tertiary
lincom secondary-tertiary

local var Mexp_binVDem_11after
reghdfe y2     ///
c.`var'##i.edu    $agg##i.edu $cov##c.`var'   ,a(country i.birthyear) cluster(country)  keepsing
local r2=e(r2_a)
eststo m3: expint y2  `var'  beta
estadd local ctrl2 "√"
estadd local ctrl1 "√"
estadd local rfe "√"
estadd local bfe "√"
estadd scalar r2a=`r2'
lincom primary-secondary
lincom primary-tertiary
lincom secondary-tertiary

local var Mexp_binVDem_6after
reghdfe y2     ///
c.`var'##i.edu    $agg##i.edu $cov##c.`var'   ,a(country i.birthyear) cluster(country)  keepsing
local r2=e(r2_a)
eststo m4: expint y2  `var'  beta
estadd local ctrl2 "√"
estadd local ctrl1 "√"
estadd local rfe "√"
estadd local bfe "√"
estadd scalar r2a=`r2'
lincom primary-secondary
lincom primary-tertiary
lincom secondary-tertiary

esttab m* using TA8.rtf , b(3) se(3) nonote nobaselevel  unstack noomitted  eqlabel(none) label  replace  star(+ 0.1 * 0.05 ** 0.01 *** 0.01) refcat(primary "{\i Effect of cumulative exposure to media freedom on}" , nolabel) /// 
varlabel(primary "Primary education or below" secondary "Secondary education" tertiary "College education or above") ///
mgroup("{\i DV = Trust in Science (IRT)}", pattern(1 0 0 0 0 0 0 0 0 0 0) ) ///
mtitles("{\i Continuous Media Freedom Exposure (V-Dem), age 11 to present}" "{\i Continuous Media Freedom Exposure (V-Dem), age 6 to present}" "{\i Binary Media Freedom Exposure (V-Dem), age 11 to present}" "{\i Binary Media Freedom Exposure (V-Dem), age 6 to present}") ///
stats(rfe bfe ctrl1 ctrl2 r2a N, labels("Country FE" "Birth year FE" "Individual-level controls × Media freedom"  "Country-level controls × Education level" "Adjusted R²" "Observations") fmt( 0 0 0 0 3 0 ) ) ///
title({\b Robustness for Fixed Effects Model: Alternative Measures of Media Freedom Exposure and Education Level}) ///
note("{\i \b Notes:} + {\i p} < 0.1, * {\i p} < 0.05, ** {\i p} < 0.01,*** {\i p} < 0.001 (two-tailed test)." )

**#Table A.9: Subgroup Results from Multilevel Model: Media Freedom Exposure
eststo clear

local var Mexp_conVDem_14after
eststo h1: mixed y2 edu female emp income loggdppc logpop log_nobel_pc log_qs500_pc e_pelifeex `var'_std || country: if v2x_regime==0|v2x_regime==1 , mle
estat icc
estadd scalar ICC = r(icc2)*100 
qui tab country if `var' != . & v2x_regime==0|v2x_regime==1
estadd scalar N2 = r(r)

local var Mexp_conVDem_14after
eststo h2: mixed y2 edu female emp income loggdppc logpop log_nobel_pc log_qs500_pc e_pelifeex `var'_std || country: if v2x_regime==3|v2x_regime==4, mle 
estat icc
estadd scalar ICC = r(icc2)*100 
qui tab country if `var' != . & v2x_regime==3|v2x_regime==4
estadd scalar N2 = r(r)

local var Mexp_conVDem_14after
eststo h3: mixed y2 edu female emp income loggdppc logpop log_nobel_pc log_qs500_pc e_pelifeex `var'_std || country: if wbi>=3 , mle 
estat icc
estadd scalar ICC = r(icc2)*100 
qui tab country if `var' != . & wbi>=3 
estadd scalar N2 = r(r)

local var Mexp_conVDem_14after
eststo h4: mixed y2 edu female emp income loggdppc logpop log_nobel_pc log_qs500_pc e_pelifeex `var'_std || country: if wbi<3 , mle 
estat icc
estadd scalar ICC = r(icc2)*100 
qui tab country if `var' != . & wbi<3 
estadd scalar N2 = r(r)

local var Mexp_conVDem_14after
eststo h5: mixed y2 edu female emp income loggdppc logpop log_nobel_pc log_qs500_pc e_pelifeex `var'_std || country: if v2x_corr >=0.5 & v2x_corr!=., mle 
estat icc
estadd scalar ICC = r(icc2)*100 
qui tab country if `var' != . & v2x_corr >=0.5 & v2x_corr!=.
estadd scalar N2 = r(r)

local var Mexp_conVDem_14after
eststo h6: mixed y2 edu female emp income loggdppc logpop log_nobel_pc log_qs500_pc e_pelifeex `var'_std || country: if v2x_corr <0.5 & v2x_corr!=., mle 
estat icc
estadd scalar ICC = r(icc2)*100 
qui tab country if `var' != . & v2x_corr <0.5 & v2x_corr!=.
estadd scalar N2 = r(r)

esttab h* using TA9.rtf , b(3) se(3) nonote nobaselevel noomitted  eqlabel(none) label  replace  star(+ 0.1 * 0.05 ** 0.01 *** 0.01)  ///
refcat(female "{\i Individual-Level Variables}" loggdppc "{\i Country-Level Variables}" c.Mexp_conVDem_14after#c.edu "{\i Moderation effect}", nolabel) /// 
keep( _cons edu female emp income edu loggdppc logpop log_nobel_pc log_qs500_pc e_pelifeex  Mexp_conVDem_14after_std ) ///
order(_cons edu)  ///
varlabel(_cons "Constant" female "Female" emp "Employment status" income "Income" edu "Education level" loggdppc "Log GDP per capita" logpop "Log population" log_nobel_pc "Log number of Nobel Prize winners per capita" log_qs500_pc "Log number of QS 500 universities per capita" e_pelifeex "Life expectancy"  Mexp_conVDem_14after_std "Media freedom exposure" ) ///
mgroup("{\i DV = Trust in Science (IRT)}") ///
mtitles("{\i Autocratic Regime}" "{\i Democratic Regime}" "{\i High Income}" "{\i Low Income}" "{\i High Political Corruption}" "{\i Low Political Corruption}" ) ///
title({\b Subgroup Analysis of Multilevel Model: Media Freedom Exposure}) ///
stats( ICC N2 N, labels( "ICC (%)" "Country-Level Observations" "Individual-Level Observations") fmt( 2 0 0 ) ) ///
note("{\i \b Notes:} + {\i p} < 0.1, * {\i p} < 0.05, ** {\i p} < 0.01,*** {\i p} < 0.001 (two-tailed test). This table displays results from several subgroup analyses based on the multilevel liner models with country-level random intercepts, with standardized media freedom exposure as the independent variable. The first two columns present results from democracies and non-democracies respectively, classified by the Regimes of the World index in 2020 from the V-Dem dataset. Columns 3 and 4 present results subset by country income levels from World Bank's definition in 2020. Columns 5 and 6 present the results subset by political corruption, specifically based on whether the 2020 Political Corruption index value from the V-Dem dataset was above or below 0.5. Specifications are based on Column 4 of Tables 1.")

**#Table A.10: Subgroup Results from Fixed Effects Model: Media Freedom Exposure and Education Level
eststo clear

local var Mexp_conVDem_14after
reghdfe y2     ///
c.`var'##i.edu    $agg##i.edu $cov##c.`var'  if v2x_regime==0|v2x_regime==1 ,a(country i.birthyear) cluster(country)  keepsing
local r2=e(r2_a)
eststo m1: expint y2  `var'  beta
estadd scalar r2a=`r2'

local var Mexp_conVDem_14after
reghdfe y2     ///
c.`var'##i.edu    $agg##i.edu $cov##c.`var'  if v2x_regime==3|v2x_regime==4 ,a(country i.birthyear) cluster(country)  keepsing
local r2=e(r2_a)
eststo m2: expint y2  `var'  beta
estadd scalar r2a=`r2'

local var Mexp_conVDem_14after
reghdfe y2     ///
c.`var'##i.edu    $agg##i.edu $cov##c.`var'  if wbi>=3 & wbi!=. ,a(country i.birthyear) cluster(country)  keepsing
local r2=e(r2_a)
eststo m3: expint y2  `var'  beta
estadd scalar r2a=`r2'

local var Mexp_conVDem_14after
reghdfe y2     ///
c.`var'##i.edu    $agg##i.edu $cov##c.`var'  if wbi<3 & wbi!=. ,a(country i.birthyear) cluster(country)  keepsing
local r2=e(r2_a)
eststo m4: expint y2  `var'  beta
estadd scalar r2a=`r2'

local var Mexp_conVDem_14after
reghdfe y2     ///
c.`var'##i.edu    $agg##i.edu $cov##c.`var'  if v2x_corr >=0.5 & v2x_corr!=. ,a(country i.birthyear) cluster(country)  keepsing
local r2=e(r2_a)
eststo m5: expint y2  `var'  beta
estadd scalar r2a=`r2'

local var Mexp_conVDem_14after
reghdfe y2     ///
c.`var'##i.edu    $agg##i.edu $cov##c.`var'  if v2x_corr <0.5 & v2x_corr!=. ,a(country i.birthyear) cluster(country)  keepsing
local r2=e(r2_a)
eststo m6: expint y2  `var'  beta
estadd scalar r2a=`r2'

estadd local ctrl2 "√", replace: m*
estadd local ctrl1 "√", replace: m*
estadd local rfe "√", replace: m*
estadd local bfe "√", replace: m*

esttab m* using TA10.rtf , b(3) se(3) nonote nobaselevel  unstack noomitted  eqlabel(none) label  replace  star(+ 0.1 * 0.05 ** 0.01 *** 0.01) refcat(primary "{\i Effect of cumulative exposure to democracy on}" , nolabel) /// 
varlabel(primary "Primary education or below" secondary "Secondary education" tertiary "College education or above") ///
mgroup("{\i DV = Trust in Science (IRT)}", pattern(1 0 0 0 0 0 0) ) ///
mtitles("{\i Autocratic Regime}" "{\i Democratic Regime}" "{\i High Income}" "{\i Low Income}" "{\i High Political Corruption}" "{\i Low Political Corruption}") ///
stats(rfe ctrl1 ctrl2 r2a N, labels("Country and Birth year FE" "Country-level controls × Education level"  "Individual-level controls × Media freedom" "Adjusted R²"  "Observations") fmt( 0 0 0 3 0 ) ) ///
title({\b Subgroup Analysis from Fixed Effects Model: Media Freedom Exposure and Education Level}) ///
note("{\i \b Notes:} + {\i p} < 0.1, * {\i p} < 0.05, ** {\i p} < 0.01,*** {\i p} < 0.001 (two-tailed test). This table presents results from several subgroup analyses based on the fixed effects model. The subgroups are the same as those in Table A.17. Specifications are based on Column 3 of Tables 1 and 2. Standard errors clustered at country level are reported in parentheses." )

**#Table A.11: Multilevel Model: Addressing the Issue of Preference Falsification: Current Countries Whose Media Freedom Score (V-Dem) Ranking Within the Top 45%
eststo clear

local var Mexp_conVDem_14after
eststo h1: mixed y2 edu female emp income loggdppc logpop log_nobel_pc log_qs500_pc e_pelifeex `var'_std || country:, mle
estat icc
estadd scalar ICC = r(icc2)*100
qui tab country if `var'_std != .
estadd scalar N2 = r(r)

local var Mexp_binVDem_14after
eststo h2: mixed y2 edu female emp income loggdppc logpop log_nobel_pc log_qs500_pc e_pelifeex `var'_std || country:, mle
estat icc
estadd scalar ICC = r(icc2)*100
qui tab country if `var'_std != .
estadd scalar N2 = r(r)

local var Mexp_DSP_14after
eststo h3: mixed y2 edu female emp income loggdppc logpop log_nobel_pc log_qs500_pc e_pelifeex `var'_std || country:, mle
estat icc
estadd scalar ICC = r(icc2)*100
qui tab country if `var'_std != .
estadd scalar N2 = r(r)

local var Mexp_GMF_14after
eststo h4: mixed y2 edu female emp income loggdppc logpop log_nobel_pc log_qs500_pc e_pelifeex `var'_std || country:, mle
estat icc
estadd scalar ICC = r(icc2)*100
qui tab country if `var'_std != .
estadd scalar N2 = r(r)

esttab h* using TA11.rtf , b(3) se(3) nonote nobaselevel noomitted  eqlabel(none) label  replace  star(+ 0.1 * 0.05 ** 0.01 *** 0.01)  /// 
refcat( edu "{\i Individual-Level Variables}" loggdppc "{\i Country-Level Variables}" , nolabel) /// 
keep( _cons $control Mexp_conVDem_14after_std Mexp_binVDem_14after_std Mexp_DSP_14after_std Mexp_GMF_14after_std ) ///
order(_cons) ///
varlabel( _cons "Constant" female "Female" emp "Employment status" income "Income" edu "Education level" loggdppc "Log GDP per capita" logpop "Log population" log_nobel_pc "Log number of Nobel Prize winners per capita" log_qs500_pc "Log number of QS 500 universities per capita" e_pelifeex "Life expectancy" Mexp_conVDem_14after_std "Media freedom exposure") ///
mgroup("{\i DV = Trust in Science (IRT)}") ///
mtitles("{\i Continuous Media Freedom Exposure (V-Dem)}" "{\i Binary Media Freedom Exposure (V-Dem)}" "{\i Digital Media Freedom Exposure (DSP)}" "{\i Media Freedom Exposure (GMF)}") ///
title({\b Multilevel Model: Addressing the Issue of Preference Falsification of Fixed Effects Model: Current Countries Whose Media Freedom (V-Dem) Ranking Within the Top 45%}) ///
stats( ICC N2 N, labels( "ICC (%)" "Country-Level Observations" "Individual-Level Observations") fmt( 2 0 0 ) ) ///
note("{\i \b Notes:} + {\i p} < 0.1, * {\i p} < 0.05, ** {\i p} < 0.01,*** {\i p} < 0.001 (two-tailed test)." )

**#Table A.12: Fixed Effects Model: Addressing the Issue of Preference Falsification: Current Countries Whose Media Freedom Score (V-Dem) Ranking Within the Top 45%
eststo clear

local var Mexp_conVDem_14after
reghdfe y2     ///
c.`var'##i.edu    $agg##i.edu $cov##c.`var' ///
if Mf_GPCAVDem_2020 >= Mf_GPCAVDem_2020_55  ///
,a(country i.birthyear) cluster(country)  keepsing
local r2=e(r2_a)
eststo m1: expint y2  `var'  beta
estadd local ctrl2 "√"
estadd local ctrl1 "√"
estadd local rfe "√"
estadd local bfe "√"
estadd scalar r2a=`r2'

local var Mexp_binVDem_14after
reghdfe y2     ///
c.`var'##i.edu    $agg##i.edu $cov##c.`var' ///
if Mf_GPCAVDem_2020 >= Mf_GPCAVDem_2020_55   ///
,a(country i.birthyear) cluster(country)  keepsing
local r2=e(r2_a)
eststo m2: expint y2  `var'  beta
estadd local ctrl2 "√"
estadd local ctrl1 "√"
estadd local rfe "√"
estadd local bfe "√"
estadd scalar r2a=`r2'

local var Mexp_DSP_14after
reghdfe y2     ///
c.`var'##i.edu    $agg##i.edu $cov##c.`var' ///
if Mf_GPCAVDem_2020 >= Mf_GPCAVDem_2020_55   ///
,a(country i.birthyear) cluster(country)  keepsing
local r2=e(r2_a)
eststo m3: expint y2  `var'  beta
estadd local ctrl2 "√"
estadd local ctrl1 "√"
estadd local rfe "√"
estadd local bfe "√"
estadd scalar r2a=`r2'

local var Mexp_GMF_14after
reghdfe y2     ///
c.`var'##i.edu    $agg##i.edu $cov##c.`var' ///
if Mf_GPCAVDem_2020 >= Mf_GPCAVDem_2020_55  ///
,a(country i.birthyear) cluster(country)  keepsing
local r2=e(r2_a)
eststo m4: expint y2  `var'  beta
estadd local ctrl2 "√"
estadd local ctrl1 "√"
estadd local rfe "√"
estadd local bfe "√"
estadd scalar r2a=`r2'

esttab m* using TA12.rtf , b(3) se(3) nonote nobaselevel  unstack noomitted  eqlabel(none) label  replace  star(+ 0.1 * 0.05 ** 0.01 *** 0.01) refcat(primary "{\i Effect of cumulative exposure to democracy on}" , nolabel) /// 
varlabel(primary "Primary education or below" secondary "Secondary education" tertiary "College education or above") ///
mgroup("{\i DV = Trust in Science (IRT)}", pattern(1 0 0 0 0 0 0) ) ///
mtitles("{\i Continuous Media Freedom Exposure (V-Dem)}" "{\i Binary Media Freedom Exposure (V-Dem)}" "{\i Digital Media Freedom Exposure (DSP)}" "{\i Media Freedom Exposure (GMF)}") ///
stats(rfe bfe ctrl1 ctrl2 r2a N, labels("Country FE" "Birth year FE" "Country-level controls × Education level"  "Individual-level controls × Media freedom" "Adjusted R²"  "Observations") fmt( 0 0 0 0 3 0 ) ) ///
title({\b Addressing the Issue of Preference Falsification of Fixed Effects Model: Current Countries Whose Media Freedom (V-Dem) Ranking Within the Top 45%}) ///
note("{\i \b Notes:} + {\i p} < 0.1, * {\i p} < 0.05, ** {\i p} < 0.01,*** {\i p} < 0.001 (two-tailed test).")

**#Table A.13: Placebo Analysis of Multilevel Model: Media Freedom Exposure
use Data_Gu_Main, clear

keep female emp income edu loggdppc logpop log_nobel_pc log_qs500_pc e_pelifeex y2 Mexp_conVDem_14after_std country

mixed y2 female emp income edu loggdppc logpop log_nobel_pc log_qs500_pc e_pelifeex Mexp_conVDem_14after_std || country:, mle
*Coefficient：-0.091756
*p-value : 0.000
*z-value : -9.43

permute Mexp_conVDem_14after_std beta = _b[Mexp_conVDem_14after_std] se = _se[Mexp_conVDem_14after_std] df = e(df_m), ///
reps(1000) rseed(123) saving("HLM_PlaceboAnalysis.dta", replace): ///
mixed y2 female emp income edu loggdppc logpop log_nobel_pc log_qs500_pc e_pelifeex Mexp_conVDem_14after_std || country:, mle

/**回归系数图**/
use "HLM_PlaceboAnalysis.dta", clear

gen z_value = beta / se
gen p_value = 2 * (1 - normal(abs(beta/se)))

/**Coefficient**/
sum beta, detail  

twoway(kdensity beta,  ///
xline(`r(mean)', lpattern(dash)  lcolor(black))  ///
title("(a)") ///
xtitle("Coefficient")  ///
ytitle("Kernel Density")  ///
saving(placebo_HLM_Coefficient1, replace)), ///
xlabel(,format(%05.3f))  ///
ylabel(,format(%2.0f))  

*xline(-0.0886132 , lpattern(solid) lcolor(black)) 

graph export "placebo_HLM_Coefficient1.png", replace

/**P值**/
sum beta, detail

twoway(scatter p_value beta,  ///
msy(oh) mcolor(black)  ///
title("(b)") ///
xline(`r(mean)' , lpattern(dash) lcolor(black)) ///
yline(0.000 , lpattern(shortdash) lcolor(black))  ///
xtitle("Coefficient") ///
ytitle("p-value")  ///
saving(placebo_HLM_Pvalue1, replace)),  ///
ylabel(0(0.25)1, format(%03.2f))

*xline(-0.0886132 , lpattern(dash) lcolor(black))

graph export "placebo_HLM_Pvalue1.png", replace

/**t值**/
sum z_value, detail

twoway(kdensity z_value, ///
xline(`r(mean)' , lpattern(dash) lcolor(black))  ///
title("(c)") ///
xtitle("z-value") ///
ytitle("Kernel Density")  ///
saving(placebo_HLM_Zvalue1, replace)),  ///
xlabel(,format(%01.0f))  ///
ylabel(,format(%02.1f)) 

*xline(-5.36 , lpattern(solid) lcolor(black))

graph export "placebo_HLM_Zvalue1.png", replace

*合并三个图
graph combine placebo_HLM_Coefficient1.gph placebo_HLM_Pvalue1.gph placebo_HLM_Zvalue1.gph, ///
xsize(10) ysize(15) cols(1)

graph export "placebo_HLM_All.png", replace

**#Table A.14: Placebo Analysis of Fixed Effects Model: Media Freedom Exposure and Education Level
use Data_Gu_Main, clear

keep female emp income edu_nc loggdppc logpop log_nobel_pc log_qs500_pc e_pelifeex y2 Mexp_conVDem_14after_std Mexp_conVDem_14after_stdXedu_nc country birthyear

global cov ( i.female i.emp i.income )
global agg ( c.loggdppc c.logpop c.log_nobel_pc c.log_qs500_pc c.e_pelifeex )

reghdfe y2  Mexp_conVDem_14after_stdXedu_nc  $cov##c.Mexp_conVDem_14after_std  $agg##i.edu_nc ,a(country i.birthyear) cluster(country)  keepsing
*Coefficient：-0.064753
*p-value : 0.000
*z-value : -4.29

permute Mexp_conVDem_14after_stdXedu_nc beta = _b[Mexp_conVDem_14after_stdXedu_nc] se = _se[Mexp_conVDem_14after_stdXedu_nc] df = e(df_m), ///
reps(1000) rseed(123) saving("FE_PlaceboAnalysis.dta", replace): ///
reghdfe y2  Mexp_conVDem_14after_stdXedu_nc  $cov##c.Mexp_conVDem_14after_std  $agg##i.edu_nc ,a(country i.birthyear) cluster(country)  keepsing

/**回归系数图**/
use "FE_PlaceboAnalysis.dta", clear

gen z_value = beta / se
gen p_value = 2 * (1 - normal(abs(beta/se)))

/**Coefficient**/
sum beta, detail  

twoway(kdensity beta,  ///
xline(`r(mean)', lpattern(dash)  lcolor(black))  ///
title("(a)") ///
xtitle("Coefficient")  ///
ytitle("Kernel Density")  ///
saving(placebo_HLM_Coefficient1, replace)), ///
xlabel(,format(%05.3f))  ///
ylabel(,format(%2.0f))  

*xline(-0.0886132 , lpattern(solid) lcolor(black)) 

graph export "placebo_FE_Coefficient1.png", replace

/**P值**/
sum beta, detail

twoway(scatter p_value beta,  ///
msy(oh) mcolor(black)  ///
title("(b)") ///
xline(`r(mean)' , lpattern(dash) lcolor(black)) ///
yline(0.000 , lpattern(shortdash) lcolor(black))  ///
xtitle("Coefficient") ///
ytitle("p-value")  ///
saving(placebo_HLM_Pvalue1, replace)),  ///
ylabel(0(0.25)1, format(%03.2f))

*xline(-0.0886132 , lpattern(dash) lcolor(black))

graph export "placebo_FE_Pvalue1.png", replace

/**t值**/
sum z_value, detail

twoway(kdensity z_value, ///
xline(`r(mean)' , lpattern(dash) lcolor(black))  ///
title("(c)") ///
xtitle("z-value") ///
ytitle("Kernel Density")  ///
saving(placebo_HLM_Zvalue1, replace)),  ///
xlabel(,format(%01.0f))  ///
ylabel(,format(%02.1f)) 

*xline(-5.36 , lpattern(solid) lcolor(black))

graph export "placebo_FE_Zvalue1.png", replace

*合并三个图
graph combine placebo_HLM_Coefficient1.gph placebo_HLM_Pvalue1.gph placebo_HLM_Zvalue1.gph, ///
xsize(10) ysize(15) cols(1)

graph export "placebo_FE_All.png", replace

**#Table A.15: Multilevel Model: Media Freedom Exposure and Science Literacy
use Data_Gu_Main, clear

eststo clear
local var Mexp_conVDem_14after
eststo h1: mixed y2 edu female emp income loggdppc logpop log_nobel_pc log_qs500_pc e_pelifeex c.`var'_std##c.science_literacy || country:, mle
estat icc
estadd scalar ICC = r(icc2)*100
qui tab country if `var' != .
estadd scalar N2 = r(r)

local var Mexp_binVDem_14after
eststo h2: mixed y2 edu female emp income loggdppc logpop log_nobel_pc log_qs500_pc e_pelifeex c.`var'_std##c.science_literacy || country:, mle
estat icc
estadd scalar ICC = r(icc2)*100
qui tab country if `var' != .
estadd scalar N2 = r(r)

local var Mexp_DSP_14after
eststo h3: mixed y2 edu female emp income loggdppc logpop log_nobel_pc log_qs500_pc e_pelifeex c.`var'_std##c.science_literacy || country:, mle
estat icc
estadd scalar ICC = r(icc2)*100
qui tab country if `var' != .
estadd scalar N2 = r(r)

local var Mexp_scoreRSF_14after 
eststo h4: mixed y2 edu female emp income loggdppc logpop log_nobel_pc log_qs500_pc e_pelifeex c.`var'_std##c.science_literacy || country:, mle
estat icc
estadd scalar ICC = r(icc2)*100
qui tab country if `var' != .
estadd scalar N2 = r(r)

local var Mexp_conFOTN_14after 
eststo h5: mixed y2 edu female emp income loggdppc logpop log_nobel_pc log_qs500_pc e_pelifeex c.`var'_std##c.science_literacy || country:, mle
estat icc
estadd scalar ICC = r(icc2)*100
qui tab country if `var' != .
estadd scalar N2 = r(r)

local var Mexp_FOTPscore_14after
eststo h6: mixed y2 edu female emp income loggdppc logpop log_nobel_pc log_qs500_pc e_pelifeex c.`var'_std##c.science_literacy || country:, mle
estat icc
estadd scalar ICC = r(icc2)*100
qui tab country if `var' != .
estadd scalar N2 = r(r)

local var Mexp_GMF_14after
eststo h7: mixed y2 edu female emp income loggdppc logpop log_nobel_pc log_qs500_pc e_pelifeex c.`var'_std##c.science_literacy || country:, mle
estat icc
estadd scalar ICC = r(icc2)*100
qui tab country if `var' != .
estadd scalar N2 = r(r)

esttab h* using TA15.rtf , b(3) se(3) nonote nobaselevel noomitted  eqlabel(none) label  replace  star(+ 0.1 * 0.05 ** 0.01 *** 0.01)  /// 
refcat(edu "{\i Individual-Level Variables}" loggdppc "{\i Country-Level Variables}" , nolabel) /// 
keep( _cons $control science_literacy Mexp_conVDem_14after_std Mexp_binVDem_14after_std Mexp_DSP_14after_std Mexp_scoreRSF_14after_std Mexp_conFOTN_14after_std Mexp_FOTPscore_14after_std Mexp_GMF_14after_std c.Mexp_conVDem_14after_std#c.science_literacy c.Mexp_binVDem_14after_std#c.science_literacy c.Mexp_DSP_14after_std#c.science_literacy c.Mexp_scoreRSF_14after_std#c.science_literacy c.Mexp_conFOTN_14after_std#c.science_literacy c.Mexp_FOTPscore_14after_std#c.science_literacy c.Mexp_GMF_14after_std#c.science_literacy ) ///
order(_cons edu female emp income science_literacy) ///
varlabel( _cons "Constant" female "Female" emp "Employment status" income "Income" edu "Education level" science_literacy "Science literacy" loggdppc "Log GDP per capita" logpop "Log population" log_nobel_pc "Log number of Nobel Prize winners per capita" log_qs500_pc "Log number of QS 500 universities per capita" e_pelifeex "Life expectancy" Mexp_conVDem_14after_std "Media freedom exposure" c.Mexp_conVDem_14after_std#c.science_literacy "Media freedom exposure × Science Literacy") ///
mgroup("{\i DV = Trust in Science (IRT)}") ///
mtitles("{\i Continuous Media Freedom Exposure (V-Dem)}" "{\i Binary Media Freedom Exposure (V-Dem)}" "{\i Digital Media Freedom Exposure (DSP)}" "{\i Press Freedom Exposure (RSF)}" "{\i Internet Freedom Exposure (FOTN)}" "{\i Press Freedom Exposure (FOTP)}" "{\i Media Freedom Exposure (GMF)}") ///
title({\b Multilevel Model: Media Freedom Exposure and Science Literacy }) ///
stats( ICC N2 N, labels( "ICC (%)" "Country-Level Observations" "Individual-Level Observations") fmt( 2 0 0 ) ) ///
note("{\i \b Notes:} + {\i p} < 0.1, * {\i p} < 0.05, ** {\i p} < 0.01,*** {\i p} < 0.001 (two-tailed test)." )

**#Table A.16: Fxied Effects Model: Media Freedom Exposure and Science Literacy
eststo clear

local var Mexp_conVDem_14after_std
local mod science_literacy
eststo n1: reghdfe y2     ///
c.`var'##c.`mod'   $cov##c.`var'    $agg##c.`mod'  ,a(country   i.birthyear##i.edu) cluster(country)  keepsing

local var Mexp_binVDem_14after_std
local mod science_literacy
eststo n2: reghdfe y2     ///
c.`var'##c.`mod'   $cov##c.`var'    $agg##c.`mod'  ,a(country   i.birthyear##i.edu) cluster(country)  keepsing

local var Mexp_DSP_14after_std
local mod science_literacy
eststo n3: reghdfe y2     ///
c.`var'##c.`mod'   $cov##c.`var'    $agg##c.`mod'  ,a(country   i.birthyear##i.edu) cluster(country)  keepsing

local var Mexp_scoreRSF_14after_std
local mod science_literacy
eststo n4: reghdfe y2     ///
c.`var'##c.`mod'   $cov##c.`var'    $agg##c.`mod'  ,a(country   i.birthyear##i.edu) cluster(country)  keepsing

local var Mexp_conFOTN_14after_std 
local mod science_literacy
eststo n5: reghdfe y2     ///
c.`var'##c.`mod'   $cov##c.`var'    $agg##c.`mod'  ,a(country   i.birthyear##i.edu) cluster(country)  keepsing

local var Mexp_FOTPscore_14after_std
local mod science_literacy
eststo n6: reghdfe y2     ///
c.`var'##c.`mod'   $cov##c.`var'    $agg##c.`mod'  ,a(country   i.birthyear##i.edu) cluster(country)  keepsing

local var Mexp_GMF_14after_std
local mod science_literacy
eststo n7: reghdfe y2     ///
c.`var'##c.`mod'   $cov##c.`var'    $agg##c.`mod'  ,a(country   i.birthyear##i.edu) cluster(country)  keepsing

estadd local ctrl1 "√", replace: n*
estadd local ctrl2 "√", replace: n*
estadd local rfe "√", replace: n*
estadd local bfe "√", replace: n*

esttab n* using TA16.rtf , b(3) se(3) nonote nobaselevel noomitted eqlabel(none) label replace star(+ 0.1 * 0.05 ** 0.01 *** 0.01)  /// 
keep (science_literacy Mexp_conVDem_14after_std Mexp_binVDem_14after_std Mexp_DSP_14after_std Mexp_scoreRSF_14after_std Mexp_conFOTN_14after_std Mexp_FOTPscore_14after_std Mexp_GMF_14after_std c.Mexp_conVDem_14after_std#c.science_literacy c.Mexp_binVDem_14after_std#c.science_literacy c.Mexp_DSP_14after_std#c.science_literacy c.Mexp_scoreRSF_14after_std#c.science_literacy c.Mexp_conFOTN_14after_std#c.science_literacy c.Mexp_FOTPscore_14after_std#c.science_literacy c.Mexp_GMF_14after_std#c.science_literacy ) ///
varlabel(science_literacy "Science Literacy" Mexp_conVDem_14after_std "Media Freedom Exposure" c.Mexp_conVDem_14after_std#c.science_literacy "Media Freedom Exposure × science literacy" ) ///
order(Mexp_conVDem_14after_std science_literacy c.Mexp_conVDem_14after_std#c.science_literacy) ///
mgroup("{\i DV = Trust in Science (IRT)}") ///
mtitles("{\i Continuous Media Freedom Exposure (V-Dem)}" "{\i Binary Media Freedom Exposure (V-Dem)}" "{\i Digital Media Freedom Exposure (DSP)}" "{\i Press Freedom Exposure (RSF)}" "{\i Internet Freedom Exposure (FOTN)}" "{\i Press Freedom Exposure (FOTP)}" "{\i Media Freedom Exposure (GMF)}") ///
stats(rfe bfe ctrl1 ctrl2 r2 N, labels("Country FE" "Birth year FE" "Country-level controls × Science Literacy"  "Individual-level controls × Media freedom" "Adjusted R²"  "Observations") fmt( 0 0 0 0 3 0 ) ) ///
title({\b Fixed Effects Model: Media Freedom Exposure and Science Literacy}) ///
note("{\i \b Notes:} + {\i p} < 0.1, * {\i p} < 0.05, ** {\i p} < 0.01,*** {\i p} < 0.001 (two-tailed test)." )

**#Table A.17: Alternative Model of Multilevel Analysis: Media Freedom Exposure and Education Level
eststo clear

local var Mexp_conVDem_14after
eststo h1: mixed y2 edu female emp income loggdppc logpop log_nobel_pc log_qs500_pc e_pelifeex `var'_std i.birthyear || country_new: 
estat icc
estadd scalar ICC = r(icc2)*100
qui tab country if `var'_std != .
estadd scalar N2 = r(r)
estadd local ri "√"
estadd local bfe "√"

local var Mexp_conVDem_14after
eststo h2: mixed y2 edu female emp income loggdppc logpop log_nobel_pc log_qs500_pc e_pelifeex c.`var'_std##c.edu i.birthyear || country_new: 
estat icc
estadd scalar ICC = r(icc2)*100
qui tab country if `var'_std != .
estadd scalar N2 = r(r)
estadd local ri "√"
estadd local bfe "√"

local var Mexp_conVDem_14after
eststo h3: mixed y2 edu female emp income loggdppc logpop log_nobel_pc log_qs500_pc e_pelifeex `var'_std || country_new: `var'_std 
estat icc
estadd scalar ICC = r(icc2)*100
qui tab country if `var'_std != .
estadd scalar N2 = r(r)
estadd local ri "√"
estadd local rs "√"

local var Mexp_conVDem_14after
eststo h4: mixed y2 edu female emp income loggdppc logpop log_nobel_pc log_qs500_pc e_pelifeex c.`var'_std##c.edu || country_new: `var'_std 
estat icc
estadd scalar ICC = r(icc2)*100
qui tab country if `var'_std != .
estadd scalar N2 = r(r)
estadd local ri "√"
estadd local rs "√"

local var Mexp_conVDem_14after
eststo h5: mixed y2 edu female emp income loggdppc logpop log_nobel_pc log_qs500_pc e_pelifeex `var'_std i.birthyear || country_new: `var'_std 
estat icc
estadd scalar ICC = r(icc2)*100
qui tab country if `var'_std != .
estadd scalar N2 = r(r)
estadd local ri "√"
estadd local rs "√"
estadd local bfe "√"

local var Mexp_conVDem_14after
eststo h6: mixed y2 edu female emp income loggdppc logpop log_nobel_pc log_qs500_pc e_pelifeex c.`var'_std##c.edu i.birthyear || country_new: `var'_std 
estat icc
estadd scalar ICC = r(icc2)*100
qui tab country if `var'_std != .
estadd scalar N2 = r(r)
estadd local ri "√"
estadd local rs "√"
estadd local bfe "√"

local var Mexp_binVDem_14after
eststo h7: mixed y2 edu female emp income loggdppc logpop log_nobel_pc log_qs500_pc e_pelifeex `var'_std i.birthyear || country_new: 
estat icc
estadd scalar ICC = r(icc2)*100
qui tab country if `var'_std != .
estadd scalar N2 = r(r)
estadd local ri "√"
estadd local bfe "√"

local var Mexp_binVDem_14after
eststo h8: mixed y2 edu female emp income loggdppc logpop log_nobel_pc log_qs500_pc e_pelifeex c.`var'_std##c.edu i.birthyear || country_new: 
estat icc
estadd scalar ICC = r(icc2)*100
qui tab country if `var'_std != .
estadd scalar N2 = r(r)
estadd local ri "√"
estadd local bfe "√"

local var Mexp_binVDem_14after
eststo h9: mixed y2 edu female emp income loggdppc logpop log_nobel_pc log_qs500_pc e_pelifeex `var'_std || country_new: `var'_std 
estat icc
estadd scalar ICC = r(icc2)*100
qui tab country if `var'_std != .
estadd scalar N2 = r(r)
estadd local ri "√"
estadd local rs "√"

local var Mexp_binVDem_14after
eststo h10: mixed y2 edu female emp income loggdppc logpop log_nobel_pc log_qs500_pc e_pelifeex c.`var'_std##c.edu || country_new: `var'_std 
estat icc
estadd scalar ICC = r(icc2)*100
qui tab country if `var'_std != .
estadd scalar N2 = r(r)
estadd local ri "√"
estadd local rs "√"

local var Mexp_binVDem_14after
eststo h11: mixed y2 edu female emp income loggdppc logpop log_nobel_pc log_qs500_pc e_pelifeex `var'_std i.birthyear || country_new: `var'_std 
estat icc
estadd scalar ICC = r(icc2)*100
qui tab country if `var'_std != .
estadd scalar N2 = r(r)
estadd local ri "√"
estadd local rs "√"
estadd local bfe "√"

local var Mexp_binVDem_14after
eststo h12: mixed y2 edu female emp income loggdppc logpop log_nobel_pc log_qs500_pc e_pelifeex c.`var'_std##c.edu i.birthyear || country_new: `var'_std 
estat icc
estadd scalar ICC = r(icc2)*100
qui tab country if `var'_std != .
estadd scalar N2 = r(r)
estadd local ri "√"
estadd local rs "√"
estadd local bfe "√"

esttab h* using TA17.rtf , b(3) se(3) nonote nobaselevel noomitted  eqlabel(none) label  replace star(+ 0.1 * 0.05 ** 0.01 *** 0.01)  /// 
refcat(edu "{\i Individual-Level Variables}" loggdppc "{\i Country-Level Variables}" , nolabel) /// 
keep( _cons $control Mexp_conVDem_14after_std Mexp_binVDem_14after_std c.Mexp_conVDem_14after_std#c.edu c.Mexp_binVDem_14after_std#c.edu ) ///
order(_cons edu) ///
varlabel( _cons "Constant" female "Female" emp "Employment status" income "Income" edu "Education level" loggdppc "Log GDP per capita" logpop "Log population" log_nobel_pc "Log number of Nobel Prize winners per capita" log_qs500_pc "Log number of QS 500 universities per capita" e_pelifeex "Life expectancy" Mexp_conVDem_14after_std "Media freedom exposure" c.Mexp_conVDem_14after_std#c.edu "Media freedom exposure × Education level" ) ///
mgroup("{\i DV = Trust in Science (IRT)}") ///
mtitles("{\i Continuous Media Freedom Exposure (V-Dem)}" "" "" "" "" "" "{\i Binary Media Freedom Exposure (V-Dem)}" "" "" "" "" "" )  ///
title({\b Alternative Model of Multilevel Analysis: Media Freedom Exposure and Education Level }) ///
stats( ri rs bfe ICC N2 N, labels( "Random intercept" "Random slope" "Birth year FE" "ICC (%)" "Country-Level Observations" "Individual-Level Observations") fmt( 0 0 0 2 0 0 ) ) ///
note("{\i \b Notes:} + {\i p} < 0.1, * {\i p} < 0.05, ** {\i p} < 0.01,*** {\i p} < 0.001 (two-tailed test)." )
 